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SUMMARY 

A control -volume based finite difference method to solve the Reynolds 
averaged Navier- stokes equations is presented. A pressure correction 
equation valid at all flow velocities and a pressure staggered grid layout 
are used In the method. Example problems presented herein include: a 
developing laminar channel flow, a developing laminar pipe flow, a 
lid-driven square cavity flow, a laminar flow through a 90-degree bent 
channel, a laminar polar cavity flow, and a turbulent supersonic flow over 
a compression ramp. A k-c turbulence model supplemented with a near-wall 
turbulence model was used to solve the turbulent flow. It is shown that the 
method yields accurate computational results even when highly skewed, 
unequally spaced, curved grids are used. It is also shown that the method 
is strongly convergent for high Reynolds number flows. 


*Work funded under Space Act Agreement C99066G. 
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Nomenclature 

coefficient for incremental u-velocity 
coefficient for incremental v-velocity 
constant coefficient for f^ equation (-0.025) 
constant coefficient for f^ equation (-0.00001) 
constant coefficient for f e equation 
turbulence model constants for e equation, (i-1,2) 
constant coefficient for eddy viscosity equation (« 
normal distance from wall 

wall damping function for eddy viscosity equation 

wall damping function for € w equation 

turbulent kinetic energy 

effective thermal conductivity (-k^U^) 

the rma 1 c onduc t i v i ty 

turbulent thermal conductivity (-Cp 

outward normal vector, (-{n x , n^}) 

pressure 

production rate of turbulent kinetic energy 
gas constant. 

Reynolds number 

turbulent Reynolds number (-k^/(i/€^)) 
temperature 

friction velocity (-%/(r w /p)) 
velocity vector («{u,v}) 
cartesian coordinates (-{x,y}) 
wall coordinate (-u r d n /t/) 
dissipation rate 
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dissipation rate of turbulent kinetic energy 

dissipation rate inside the near-wall layer 

von Karman constant (-0.41) 

molecular viscosity 

effective viscosity 

turbulent viscosity 

kinematic viscosity of fluid 

turbulent eddy viscosity 

curvilinear coordinates 

density 

turbulent Prandtl number for k-equation 

turbulent Prandtl number for energy equation 

turbulent Prandtl number for e -equation 

wall shearing stress 

dissipation function for energy equation 


Superscripts 

a non-dimensional value normalized by the free stream value 

n iteration level 

* current value 

' incremental (or corrective) value 


Subscripts 

nb neighboring grid points, (-{E, W, S, N}) 

P grid point 


Mathematical symbol 

£ summation 
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INTRODUCTION 


A control -volume based finite difference method to solve the Reynolds 
averaged Navier- Stokes equations for all flow velocities is presented. The 
method is an extension of the pressure correction method (SIMPLE) which is 
used primarily to solve incompressible flows [1,2]. Numerical methods based 
on the pressure correction method have been used extensively to solve 
complex turbulent flows [3], including chemically reacting turbulent flows 

[4] , due to their strongly convergent nature. In the present study, a 
pressure correction equation which is valid at all flow velocities is used 
for numerical calculations of incompressible and compressible flows. 

Many finite difference methods to solve the compressible flow 
equations are based on the flux- splitting method. The Beam-Warming method 

[5] and the McCormack method [6] are the representatives of the 
flux- splitting methods. The flux- splitting methods were originally 
developed to solve the Euler equations and then extended to include the 
viscous term to solve the Navier- Stokes equations. The most distinguishing 
practical difference between the pressure correction methods and the 
flux- splitting methods lies in the way the diffusion term is treated. In 
the former class of methods, the diffusion term has been incorporated into 
the stiffness matrix while, in the latter class of methods, the diffusion 
term has been incorporated into the system of equations as the load vector 
term. For turbulent flows with extensive recirculation zones, the pressure 
correction methods may be numerically more stable conceptually. However, 
the pressure correction methods have mostly been used for Incompressible 
flows and the flux splitting methods have mostly been used for compressible 
flows. Therefore, definitive advantages and disadvantages of these two 
classes of methods can not be discussed with confidence as yet. 
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The original numerical method based on the SIMPLE algorithm [2] is 
used to solve the Navier- Stokes equations whose domain can be discretized 
using orthogonal grids. A fully staggered grid layout [1] has been used in 
the method. However, in many flow problems of practical Importance, the 
boundary geometries are complex and arbitrary shaped blockages may exist 
inside the flow path. A number of papers to extend the pressure correction 
methods for flows with arbitrary geometries and for compressible flows have 
appeared in recent years [7-13]. 

A number of grid layouts have been proposed and tested to Identify the 
most suitable one to solve the Navier-Stokes equations defined on 
arbitrary, complex geometries. In Reference 7, the standard fully staggered 
grid layout, Figure l-(a), was used to solve the Navier-Stokes equations 
defined on curved geometries. This grid layout can not be used to solve 
flows inside 90-degree bent ducts (see Reference 8 for details) . A 
collocated grid layout, Figure l-(b), was used in References 9 and 10. In 
Reference 9, the velocity - pressure decoupling was prevented by including 
an artificial viscosity, while In Reference 10, the same purpose was 
achieved by evaluating the incremental velocities at mid- sides of the 
control volume. In Reference 8, the standard fully staggered grid was used 
and the velocity vector was located at all grid points except at the 
pressure grid point (see Figure l-(c)). In this case, the number of degrees 
of freedom for velocity is doubled while that of pressure remains the same 
as in the original case. Note that the accuracy of numerical solutions 
depends not only on the number of velocity grid points but also on the 
number of pressure grid points. Hence the accuracy can not be improved as 
much as the doubled number of velocity grid points might suggest. In 
References 12 and 13, the velocities were located at the same grid points 


5 



and the pressure was located at the centroid of the cell formed by the four 
adjacent velocity grid points (see Figure l-(d)). This grid layout has been 
used successfully in penalty finite element methods for a long time [14]. 

It was first used in the control -volume based finite difference method in 
Vanka et. al . [12]. They mentioned that it was not easy to obtain 

convergent solutions due to the velocity-pressure decoupling. The mechanism 
that yields the velocity-pressure decoupled solution was heuris tically 
shown in Reference [8], In Reference 13, the velocity-pressure decoupling 
was eliminated by using a non- conforming domain for the mass imbalance 
calculation. In the present study, the velocity-pressure decoupling was 
eliminated by moving the off-diagonal terms to the load vector term. The 
resulting system of equations was solved using the Tri-Diagonal Matrix 
Algorithm (TDMA) [1] . Thus any uncertainty that may arise due to the use of 
a non- conforming domain for the mass imbalance calculation does not exist 
in the present method. 

A few different numerical procedures have also been used to solve the 
pressure correction equation for compressible flows. The SIMPLE-R [1] and 
the SIMPLE-C [15] were used In References 11 and 13, respectively. The 
pressure, velocity, and density were corrected based on the Incremental 
pressure (or pressure correction) in References 11 and 13. In the present 
study, only the pressure and velocity were corrected from the incremental 
pressure as In the standard SIMPLE method [1]. Density was obtained from 
the equation of state for perfect gas so that the same numerical procedure 
could equally be applicable for numerical computations of chemically 
reacting turbulent flows In the future. The numerical procedure to solve 
the pressure correction equation for compressible flows may need to be 
studied further In the future. 
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A k-e turbulence model supplemented with a near-wall turbulence model 
was used in the present study. Establishment of the near-wall turbulence 
model and its application to fully developed turbulent channel and pipe 
flows can be found in Reference 16. It has been shown in the reference that 
the near-wall turbulence model can resolve the over-shoot phenomena of the 
turbulent kinetic energy and the dissipation rate in the region very close 
to the wall. Thus significantly improved computational results for the 
near-wall turbulence structure were obtained. It is also shown in Reference 
17 that the turbulence model yields correct location of the shock for a 
transonic flow over an axisymmetric curved hill with shock wave - turbulent 
boundary layer interaction [18]. 

A number of flow cases have been solved to test the accuracy and the 
convergence nature of the present numerical method. The example flows 
presented herein include: a developing laminar channel flow, a developing 
laminar pipe flow, a lid-driven square cavity flow [19-21], a 
two-dimensional laminar flow in a 90-degree bent channel, a polar cavity 
flow [22], and a turbulent supersonic flow over a compression ramp [23-24]. 

REYNOLDS AVERAGED NAVI ER- STOKES EQUATIONS FOR COMPRESSIBLE FLOWS 
The compressible turbulent flow equations are given as, 

d d ... 

— ( pn ) + — (pv) - 0. 
dx dy 


d d d d dp 

— (puu) + — (/>vu) ( r xx) + — ( r xy) 

dx dy dx dy dx 
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and the density is obtained from the perfect gas law given as p-/?RT. The 
turbulent Prandtl number (a T ) of 0.75 was used for the energy equation in 
the present calculations. The compressible turbulent flow equations for 
axisymmetric case can be found in Reference 17. 

The molecular viscosity and the thermal conductivity were obtained 
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from the Sutherland's laws given as [25]; 


T ' 
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where /i Q - 1.716 x 10“-* Kg/m-sec, T c - 273.1° K, S - 110.6° K; and 
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where k Q - 0.0264 Kg/m-K, T Q - 273.1° K, and S - 194.4° K. 


TURBULENCE EQUATIONS 

A k-e turbulence model supplemented with a near-wall turbulence model 
is described below [16,17], The turbulent kinetic energy equation for the 
entire flow domain is given as ; 
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where the production rate of turbulent kinetic energy (Pr) is the same as 
the dissipation function for the energy equation ($) . 

The dissipation rate inside the near-wall layer is given as; 
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The dissipation rate given as eq. (8) is used for eq. (7) in the 

near-wall region. The dissipation rate for the rest of the flow domain is 

obtained by solving the convection-diffusion equation for the dissipation 

0 

rate equation given as; 
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The turbulence model constants used are given as: 0^-0. 75, er £ -1.15, 
Cl _1.39 t and C2“1.88. These turbulence model constants approximately 
satisfy the near-wall equilibrium turbulence condition and an 
experimentally observed decay rate of the grid turbulence [26], Further 
discussion on the establishment of these constants can be found in 
References 27-28. 

The eddy viscosity inside the near-wall layer is given as; 
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( 11 ) 


u t ~ c /lf ^/i 


where f^-1 -exp(A^yR t + A 2 R t 2 ) * The wal1 damping function is a linear 
function of the distance from the wall in the viscous sublayer and becomes 
unity in the fully turbulent region. The eddy viscosity, given as eq. (11), 
grows in proportion to the cubic power of the distance from the wall. The 
eddy viscosity in the rest of the flow domain is given as; 



€ 


The partition between the near-wall region and the fully turbulent outer 
region can be located between y + greater than 100 and less than 300 
approximately. Details on the k-e turbulence model supplemented with the 
near-wall turbulence model can be found in References 16 and 17. 


NUMERICAL METHOD 

In the present method, all flow variables, except pressure, have been 
located at the same grid points and the pressure node has been located at 
the centroid of the cell formed by the four neighboring velocity grid 
points. Note that In the control -volume methods based on pressure 
correction methods, the discrete system of equations is derived by 
integrating the governing differential equations over the control volume 
[1] . For flows with arbitrary geometries, the number of interpolations to 
obtain flow variables at the cell boundaries for the present grid layout Is 
as small as for any of the grid layouts discussed previously. Enhanced 
convergence rate is partly attributed to the grid layout which requires 
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fewer interpolations . 

The pressure correction equation for compressible flows is described 
below. As in the standard pressure correction method [1], the density, the 
velocities, and the pressure are decomposed as; 


p - p + p 


TK , * 

u - u + u , 


V - V + V , 


X , . 

p - p + p 


( 13 ) 

( 14 ) 

( 15 ) 

( 16 ) 


where the superscript * denotes the current values of the flow variables 
which may not satisfy the conservation of mass equation yet. The discrete 
momentum equation for u-velocity can be written as; 

dp 

ApUp - X A^u^ + Sv u (17) 

dx 

where Ap is the coefficient of the u-velocity at the grid point P, Sv u is 
the load vector originating from the curvilinear grid structure, and the 
pressure gradient was left in continuous form deliberately. Substituting 
eqs. (13-16) into eqs. (17) yields; 

3(p*+p') 

A p (u* p +u'p) - X An b (u* nb +u' nb ) + Sv u (18) 

dx 

The discrete u-momentum equation based on the current flow variables which 
may not satisfy the conservation of mass equation can be written as; 
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( 19 ) 


dp* 

A p u *p " I A nb u nb + Sv u 

dx 

Subtracting eq. (19) from eq. (18) yields; 

dp' 

u' « - A u — ( 2 °) 

dx 

where A^— 1/Ap. In deriving eq. (20), the summation over the neighboring 
grid points and the load vectors in eqs . (18) and (19) have been neglected. 
The relationship between the incremental v-velocity and the pressure 
gradient in the y-coordinate direction can be obtained by the same 
procedure and is given as ; 

dp' 

v' - - (21) 

ay 

The Incremental pressure is related to the Incremental density as; 

p' - p ’ RT (22) 

where eq. (22) has been obtained from the equation of state. The 
conservation of mass equation is given as; 

V* (^V) - 0 


or, 


v-(p'V*) + v(p*V') + v-(p'V') - -v- (p*V*) 


( 23 ) 
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Substituting eqs . (20-22) into (23) yields, after some rearrangement; 
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( 24 ) 


where the higher order perturbation terra V(p'V') has been neglected in 
deriving eq. (24). The last term in eq. (24) represents the mass imbalance. 
Integrating eq. (24) over a pressure control volume yields; 
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T^+V Tly 

RT 
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3x dy d£ 



■ r di) dn apO 

O <{(p*A u — n* + p*Ay— iO — >ds - /-V-(p*V*)dx (25) 

, [ ax dy drj ) 


where the Green- Gauss theorem has been made use of to invert the volume 
integration into a surface integration. The discrete control volume 
equation for eq. (25) can be written as; 

App'p - A E p' E + Ayp' w + A s p' s + A N p' N + Sc p . + Sv p . (26) 

where Sc p » is a load vector of the mass imbalance and 
SVp ,-A^3p' /dft-A^p'/dq contains all the contributions made by the 
curvilinear grid. Note that the variable load vector Sv p « is a null vector 
in the first sweep. After the first sweep, the Sv p » term was updated in 
each sweep using the incremental pressure obtained in the previous sweep. 

The other flow equations were solved by the same procedure as that of 
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the pressure correction equation. However, the load vector Sv^ 

(i“U,v, t,k,0 is not a null vector and was evaluated only once in each 
iteration. The upwind difference approximation [1] has been used to solve 
the pressure correction equation; and the power law difference 
approximation [1], for the other flow equations. The incremental pressure 
is obtained by solving eq. (26), and the corresponding incremental 
velocities are obtained from eqs . (20-21). The flow variables are updated 
by using eqs. (14-16), and these updated flow variables are used in 
computing the new current flow variables by solving eqs. (2-4) together 
with the turbulence equations. The discrete finite difference system of 
equations was solved iteratively using the TDMA until the residuals became 
smaller than the prescribed convergence criteria. Each iteration consisted 
of 7 sweeps of the pressure correction equation and 3 sweeps for the rest 
of the flow equations in the flow direction and in the transverse 
direction, respectively. The convergence criteria used are; 


R 1 - lV-pV|c < e L (27) 

C~1 

r 2 - * a n )/A" +1 | < e 2 , j-l,N, (28) 

i » j i » J i 


where N c is the number of control volumes for the pressure correction 
equation; the subscript i-{u, v, p, T, k, e) denotes each flow variable ; 
the subscript j denotes each grid point; N denotes the number of degrees of 
freedom for each flow variable; and denotes the maximum magnitude of the 
i-th flow variable. The iteration was terminated when either eq. (27) or 
eq. (28) was satisfied. 

Each flow variable was updated using an under-relaxation factor [1]. 
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The under-relaxation procedure was incorporated into the systems of 
equations for all flow equations except that of the pressure correction 
equation. The pressure correction equation was solved without 
under-relaxation; however, the pressure was obtained by adding the 
incremental pressure multiplied by an under -relaxation factor. No symptom 
that might lead to a velocity-pressure decoupled solution was observed in 
the present method. For incompressible flows, R-l.OxlO^ 5 has been used in 
eq. (24) or in eq. (25). Thus, elimination of the velocity - pressure 
decoupling mechanism is achieved by the numerical method itself and has 
nothing to do with the inclusion of the convection term into the pressure 
correction equation. 


COMPUTATIONAL RESULTS 

The numerical method described in the previous sections was tested and 
evaluated by solving a number of flow cases. Example flows presented herein 
Include: a developing laminar channel flow, a developing laminar pipe flow, 
a lid-driven square cavity flow [19-21], a two-dimensional laminar flow 
through a 90-degree bent channel, a laminar polar cavity flow [22], and a 
supersonic flow over a compression ramp with shock wave - turbulent 
boundary layer Interaction [23-24]. It is shown from the first three flow 
cases that the present method does not yield any velocity-pressure 
decoupled solution, and in fact, there was no symptom of the 
velocity-pressure decoupling for any of the flow cases considered. The rest 
of the flow cases were considered to test the accuracy, the convergence 
nature, and the applicability of the present method to flow problems of 
practical Importance. Further application of the numerical method as well 
as the turbulence model to a shock wave - turbulent boundary layer 
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Interaction in a transonic flow over an ax i symmetric curved hill [18] can 
be found in Reference 17. 

Developing Laminar Channel Flow 

A developing laminar channel flow at Re*=25 Is considered below. The 
Reynolds number is based on the inlet velocity and the channel width. The 
exit boundary was located at seven channel widths downstream of the Inlet. 
Uniform velocity was prescribed at the inlet boundary; and vanishing 
gradient boundary condition for velocities, at the exit boundary. The flow 
domain was discretized by an equally spaced mesh as well as by an unequally 
spaced mesh with 41 grid points In the flow direction and 26 In the 
transverse direction in each case. The unequally spaced mesh Is shown in 
Figure 2-(a). The convergence criteria used were e^-^-l . 0x10 ^ and the 
converged solution was obtained after approximately 410 iterations . The 
residuals at the time of convergence were R 1 -9.4xl0“ 5 and R2~l.lxl0 The 
calculated velocity profile at the exit boundary is compared with the exact 
solution in Figure 2-(b). It can be seen in the figure that the present 
method almost yields the exact solution. 

Developing Laminar Pipe Flow 

A developing laminar pipe flow at Re—25 was solved to test the 
possible existence of the velocity-pressure decoupling for the axisymmetric 
flow case. The Reynolds number is based on the Inlet velocity and the 
radius of the pipe. The exit boundary was located at seven radii downstream 
of the Inlet. An equally spaced mesh and an unequally spaced mesh with 
41x26 grid points were used as in the previous case. The unequally spaced 
mesh is shown in Figure 3-(a). The boundary conditions, the initial guess, 


17 



and the convergence criteria used are the same as in the previous channel 
flow case. The converged solution was obtained after approximately 170 
iterations and the residuals at the time of convergence were R^-4.2xl0"^ 
and R2-9.1 x 10"^. Again, it can be seen in Figure 3-(b) that the present 
method almost yields the exact solution. 

Lid-Driven Square Cavity Flow 

A lid-driven cavity flow at Re-1000 is considered below to further 
confirm that the present numerical method is free of the velocity -pressure 
decoupling mechanism. As a remark, the computational results obtained by 
the finite difference methods using fine grids can be found in References 
19 and 20 and those obtained by a finite element method can be found in 
Reference 21. The computational domain was discretized by an unequally 
spaced 81 X 81 mesh with concentration of the grid points in the near wall 
region. The boundary conditions and the initial guess used are the same as 
in Reference 21. The converged solution was obtained after approximately 
590 iterations for e^-^-l .0x10'^. The residuals at the time of convergence 
were R^-2.0xl0‘ Zf and R2“1.0xl0~^. The calculated streamline contour is 
shown in Figure 4. It can be seen in the figure that the secondary vortices 
In the bottom corners were accurately resolved by the present method. 

Two-Dimensional Laminar Flow Through a 90-Degree Bent Channel 

A two-dimensional laminar flow through a 90-degree bent channel for 
Re— 1000 Is considered below. The Reynolds number is based on the channel 
width and the bulk velocity. It is shown from this flow case that the 
present numerical method can solve flows with arbitrary geometries as 
easily as flows with rectangular geometries. Note that the numerical 
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methods adopting the fully staggered grid layout and solving for the 
cartesian velocities can not be used to solve this flow case. The inlet 
boundary was located at 5 channel widths upstream of the curved section; 
and the exit boundary, at 15 channel widths downstream of the curved 
section. The flow domain was discretized by 71 grid points in the flow 
direction and 25 in the transverse direction. The velocity profile of a 
fully developed channel flow was prescribed at the inlet boundary. The 
vanishing gradient boundary condition was used for velocities at the exit 
boundary. The convergence criteria used were e l“ e 2— 1.0x10 ^ . The converged 
solution was obtained after approximately 500 iterations. The residuals at 
the time of convergence were R^-8 . 3x10'^ and R^-l.OxlO"^. The mass flow 
rate, obtained from the prescribed inlet velocity profile, through the 
inlet boundary was 0.814419 Kg/m- sec and the calculated mass flow rate 
leaving the exit boundary was 0.814401 Kg/m- sec. Hence the relative mass 
imbalance was 2.2xl(T 5 . The grids, the calculated streamline contour, and 
the pressure contour in the vicinity of the curved section for Re— 1000 is 
shown in Figure 5. 

Laminar Polar Cavity Flow 

A laminar polar cavity flow at Re— 60 and 350 is considered below to 
test the accuracy of the present numerical method. The Reynolds number is 
based on the azimuthal velocity and the radius of curvature of the lid. The 
polar cavity is schematically shown in Figure 6. The experimental data can 
be found in Reference 22. The flow domain was discretized by 81x81 grid 
points, as in the Reference 22, with concentration of grid points in the 
corner regions. The Dirichlet boundary condition for velocities was 
prescribed at all boundaries. The convergence criteria used were 
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e ^-e2“4.0xl0' 5 . For Re-60, the converged solution was obtained after 
approximately 810 iterations and the residuals were Rj-l.SxlO -3 and 
1^2—3.9x10*^. For Re-350, the converged solution was obtained after 
approximately 820 iterations and the residuals at the time of convergence 
were Rj-7.1xlO' 3 and R 2 -3.7xlO" 5 . In each case, the required computational 
time was approximately 8 minutes for the CRAY/XMP at the NASA/LeRC. 

The calculated azimuthal and radial velocity profiles at three 
azimuthal locations for Re-60 are compared with experimental data as well 
as the computational results of Reference 22 in Figures 7-<a) and 7-(b), 
respectively. It can be seen that the present computational results and 
those obtained using the standard SIMPLE method [22] compare favorably with 
experimental data. The present method yielded a slightly better radial 
velocity profile at 9-20 degrees as shown in Figure 7 -(b). The calculated 
streamline, pressure, and the vorticity contours are shown in Figure 8. It 
can be seen from the pressure contour and the vorticity contour that the 
potential core has not been well established at Re-60. 

The calculated azimuthal and radial velocity profiles at the same 
three azimuthal locations for Re-350 are compared with experimental data as 
well as those of Reference 22 in Figure 9. It can be seen in the figure 
that both computational results are in good agreement with the experimental 
data. It is mentioned in Reference 22 that the first order differencing 
method yielded inaccurate computational results for the polar cavity flow 
at Re-350. However, it can be seen in the figure that the present method 
yielded a_s accurate computational results as those obtained using the 

me ■ ..... ... ... 

second order differencing method of Reference 22. It can be seen from the 
pressure and the vorticity contours that the potential core has been well 
established at Re-350. 
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Turbulent Supersonic Flow over a Compression Ramp 


A turbulent supersonic flow over a 24-degree compression ramp is 
considered below. The experimental data can be found in References [23-24]. 
The free stream Mach number was 2.85, the boundary layer thickness of the 
approaching supersonic flow was 0.0211 meters, and the Reynolds number 
based on the free stream condition and the boundary layer thickness was 
1 . 13xl0 6 . 

In the numerical calculation, the inlet boundary was located at 2.17 
boundary layer thicknesses upstream of the corner; and the exit boundary, 
at five boundary layer thicknesses downstream of the corner. The top 
boundary was located at seven boundary layer thicknesses away from the 
wall. The flow domain was discretized by 97 grid points in the flow 
direction and 56 in the transverse direction. The partition between the 
near-wall layer and the external region was located at approximately 4.5 
per cent of the boundary layer thickness away from the wall and 14 grid 
points were allocated Inside the near-wall layer. The grid size In the 
normal direction was increased by a factor of approximately 1.2. The inlet 
boundary condition for the tangential velocity and the turbulent kinetic 
energy were obtained from experimental data for a fully developed flat 
plate flow [29]. The non-dimensional velocity and the turbulent kinetic 
energy profiles were scaled to yield a boundary layer thickness of 0.0211 
meters at the inlet boundary. Uniform static pressure and uniform enthalpy 
were also prescribed at the inlet boundary. The no-slip boundary condition 
for velocities, vanishing turbulent kinetic energy, and a constant 
temperature which corresponds to the free stream stagnation temperature 
were prescribed at the solid wall boundary. The free stream flow condition 
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was prescribed at the top boundary, and the vanishing gradient boundary 
condition was used for all flow variables at the exit boundary. The initial 
guess was obtained by extending the inlet boundary condition in the flow 
direction. The converged solution was obtained after approximately 1400 
iterations for e^- 62 - 4 . 0x10'* . At the time of convergence and R 2 were 
3.5x10'* and 4.0x10'*, respectively. The mass flow rate through the inlet 
boundary, obtained from the prescribed inlet boundary conditions, was 
68.434 Kg/m- sec and the calculated mass flow rate leaving the exit boundary 
was 68.411 Kg/m-sec. Hence the relative mass imbalance was 3.4x10'*. The 
required computational time was approximately 18 minutes for the CRAY/XMP 
at the NASA/LeRC. 

The calculated static pressure on the wall is compared with 
experimental data as well as the computational result obtained using a 
relaxation turbulence model [24] in Figure 11. In Reference 24, several 
sets of computational results obtained using various turbulence models were 
presented. The wall pressure obtained using a relaxation turbulence model 
[24] compared most favorably with the experimental data. It can be seen In 
the figure that the present turbulence model yielded slightly compressed 
pressure distribution. The level of agreement between the experimental data 
and all the other computational results of Reference 24 was almost the same 
as that of the present computational result. 

The mean velocity profiles at s/5— -2.17, 0.0, and 2.89 are compared 
with experimental data as well as with those obtained using the relaxation 
turbulence model [24] in Figure 12, where the distance (s) has been 
measured from the corner along the surface and 5 is the boundary layer 
thickness. It can be seen that the present computational results compare 
more favorably with the experimental data than does the other computational 
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result [24]. The level of agreement between the best computational result 
in Reference 24 and the experimental data was almost the same as that of 
the present case. Note that the velocity profiles obtained using the 
relaxation model compared less favorably with the experimental data than 
those obtained using the other turbulence models [24], 

The calculated streamline contour is shown in Figure 13- (a). The 
measured flow separation zone extended from s/6~-1.44 to s/6»0.5. The 
present method yielded the flow recirculation zone extending from s/5--0.72 
to s/5-0.68. The levels of agreement between the measured flow 
recirculation zone and all the computational results, including the present 
result, were almost the same. However, the relaxation model which yielded 
the best wall pressure yielded the worst flow recirculation zone. The 
calculated static pressure contour lines are shown in Figure 13- (b), where 
the pressure has been normalized by the inlet total pressure and the 
incremental pressure between the contour lines is 0.2. The calculated 
iso-Mach lines and the turbulent kinetic energy contours are shown in 
Figures 13-(c) and 13-(d), respectively. The incremental Mach number 
between the iso-Mach lines is 0.2 in Figure 13 -(c). It has been shown in 
this example that the present computational result compared as favorably 
with the experimental data as any other computational results [24], 


CONCLUSIONS 

A control -volume based finite difference method to solve the Reynolds 
averaged Navier-Stokes equations for all flow velocities has been 
presented. 

It has been shown from the developing channel flow, the developing 
pipe flow, and the lid-driven square cavity flow that the present numerical 
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method is free of the velocity-pressure decoupling. For the channel and the 
pipe flows, the method almost yields the exact solutions. For the polar 
cavity flow, the present method yielded as accurate computational results 
as the second order differencing method [22], The turbulent supersonic flow 
over the 24-degree compression ramp [23-24] was solved using a k-e 
turbulence model supplemented with a near-wall turbulence model, In the 
method, the dissipation rate inside the near-wall region was obtained from 
an algebraic equation and that for the rest of the flow domain was obtained 
by solving the differential equation for the dissipation rate. This 
approach was found to be more advantageous than the low Reynolds number 
turbulence models since the stiff dissipation rate equation in the 
near-wall region need not be solved numerically. The computational results 
for the supersonic compression corner flow compared as favorably with the 
experimental data as any other computational results [24]. 

It has also been shown that the present numerical method yields 
accurate computational results even when highly skewed, unequally spaced, 
curved grids were used. Equally Importantly, the present method was found 
to be strongly convergent for high Reynolds number flows as well as for 
flows with complex geometries. This strongly convergent nature is 
attributed, in part, to the use of the pressure staggered grid layout. 
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(d) PRESSURE STAGGERED GRID. 
FIGURE 1. - GRID LAYOUTS. 
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(b) VELOCITY PROFILE. 

FIGURE 3, - DEVELOPING LAMINAR PIPE FLOW. 
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FIGURE 4. - STREAMLINE CONTOUR FOR LID-DRIVEN CAVITY FLOW AT Re = 1000. 
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(b) PRESSURE CONTOUR. 



(C) VORTICITY CONTOUR. 

FIGURE 8. - POLAR CAVITY FLOW FOR Re = 60. 
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(b) PRESSURE CONTOUR. 



(C) VORTICITY CONTOUR. 

FIGURE 10. - POLAR CAVITY FLOW FOR Re = 350. 
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